#Code to reproduce power analyses reported in Online Appendix 

library(tidyverse)
library(tidymodels)
library(ggpubr)

#Write function - thanks R2!
power_function<-function(alpha=.05, effect.size=.2, p = .5, N, r_pre=0) {

cval <- qnorm(1 - alpha/2) #assuming two-tailed test

if (r_pre==0)  {
the.power<-pnorm(-cval+effect.size*sqrt(p*(1-p)*N))+
  pnorm(-cval - effect.size*sqrt(p*(1-p)*N))
}
if (r_pre!=0) {
the.power <- pnorm(-cval + effect.size*sqrt(p*(1-p)*N/(1-r_pre^2))) + 
  pnorm(-cval - effect.size*sqrt(p*(1-p)*N/(1-r_pre^2)))
}  
return(list(power=the.power, samplesize=N, pretreat=r_pre))
}


############
#Sample Size Variation
############

###Post only model 

effect=.2
min_n<-50 
max_n<-1000

###Run Loop - Post-Only model 
sample<-(seq(from=min_n, to = max_n, by=1))
datalist = list()  
counter=0

mat1<-matrix(data=NA, nrow=length(sample), ncol=3)

for (i in sample) {
  counter<-counter+1
  power_calc<-power_function(N = i, effect.size=effect) 
    
  mat1[counter,1] <- i 
  mat1[counter,2] <- 0
  mat1[counter,3] <- power_calc$power
}  
mat1

final<-as_tibble(mat1)  
final<-final %>% rename(N=V1,r_pre=V2,Power=V3)
final_postonly<-final

###
#Pre-Post Values 
###

#Run Loop
sample<-(seq(from=min_n, to = max_n, by=1))
datalist = list()  
counta=0
countb=0

r_pre_input<-c(.6, .75, .9)
length_r_pre_input<-length(r_pre_input)


for(h in 1:length(r_pre_input)){
  mat1<-matrix(data=NA, nrow=length(sample), ncol=3)
  
for (i in sample) {
  counta<-counta+1
  power_calc<-power_function(N = i,effect.size=effect, r_pre=r_pre_input[h]) 
  
  mat1[counta,1] <- i 
  mat1[counta,2] <- h
  mat1[counta,3] <- power_calc$power
}  
  counta<-0 
  final<-as_tibble(mat1)  
  final<-final %>% rename(N=V1,r_pre=V2, Power=V3)
  countb<-countb+1
  datalist[[countb]]<-assign(paste0("df",h), final)

}

df_final<-bind_rows(datalist)
df_final

##Combine Post only data frame with pre-post data frame

df_final_all<-bind_rows(df_final, final_postonly)

df_final_all$r_pre<-factor(df_final_all$r_pre, 
                           levels=c(0,1,2,3),
                           labels=c("Post Only","Pre-Post (r = .6)","Pre-Post (r = .75)","Pre-Post (r = .9)"))
summary(df_final_all)

x<-ggplot(df_final_all, aes(x=N, y=Power, group=r_pre))
plot1<-x+geom_segment(x=0, xend=1000, y=.8, yend=.8)+
  geom_line(aes(color=r_pre), lwd=2.5)+
  #scale_color_manual(values=c("black","gray25","gray55","gray70"))+
  #scale_linetype_manual(values=c("solid","twodash", "longdash","dotted"))+
  #geom_line(aes(color=r_pre_cat), lwd=.5, color=c("white"), linetype=2)+
  #scale_color_manual(values=c("white","white","white","white"))+
  ylab("")+xlab("")+labs(title="Cohen's D = 0.2", size=18)+
  scale_x_continuous(breaks=seq(50, 1000, 50),expand = c(.01,.01))+
  scale_y_continuous(breaks=seq(0, 1, .10))+
  theme_bw()+
  theme(plot.title = element_text(size=18))+
  theme(legend.title = element_blank())+
  theme(legend.position="top")+
  theme(legend.text=element_text(size=17))+
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15,face="bold"))
  
plot1

plot_final<-ggarrange(plot1, plot2, 
                      ncol = 2, nrow = 1, common.legend = T, legend="top")


annotate_figure(plot_final,left = text_grob("Power", color = "Black", rot = 90, size=18),
                bottom=text_grob("Total Sample Size", size=18))


ggsave(filename="Plot - Combined.png",
       device="png",
       width=18, 
       height=10, 
       units=c("in"))


